# visualization.jl
# Comprehensive visualization tools for NSEMSolver results

using LinearAlgebra
using Printf

# Check for visualization backend availability
const HAS_PLOTS, PLOTS_MODULE = begin
    local has_plots = false
    local plots_mod = nothing
    try
        plots_mod = Base.require(Main, :Plots)
        has_plots = true
    catch
        has_plots = false
    end
    (has_plots, plots_mod)
end

const HAS_PLOTLYJS, PLOTLYJS_MODULE = begin
    local has_plotlyjs = false
    local plotlyjs_mod = nothing
    try
        plotlyjs_mod = Base.require(Main, :PlotlyJS)
        has_plotlyjs = true
    catch
        has_plotlyjs = false
    end
    (has_plotlyjs, plotlyjs_mod)
end

const HAS_MAKIE, MAKIE_MODULE = begin
    local has_makie = false
    local makie_mod = nothing
    try
        makie_mod = Base.require(Main, :Makie)
        has_makie = true
    catch
        has_makie = false
    end
    (has_makie, makie_mod)
end

const HAS_WRITEVTK, WRITEVTK_MODULE = begin
    local has_vtk = false
    local vtk_mod = nothing
    try
        vtk_mod = Base.require(Main, :WriteVTK)
        has_vtk = true
    catch
        has_vtk = false
    end
    (has_vtk, vtk_mod)
end

# Utility functions for data processing
"""
    compute_velocity_magnitude(result::NSResult) -> Array

Compute the velocity magnitude field.
"""
function compute_velocity_magnitude(result::NSResult)
    if result.w !== nothing
        # 3D case
        return sqrt.(result.u.^2 .+ result.v.^2 .+ result.w.^2)
    else
        # 2D case
        return sqrt.(result.u.^2 .+ result.v.^2)
    end
end

"""
    compute_vorticity(result::NSResult) -> Array

Compute vorticity field for 2D cases (∂v/∂x - ∂u/∂y).
"""
function compute_vorticity(result::NSResult)
    if result.w !== nothing
        @warn "3D vorticity calculation not implemented yet"
        return zeros(size(result.u))
    end
    
    # Simple finite difference approximation for vorticity
    # This is a placeholder - should be replaced with proper spectral derivatives
    nx, ny = size(result.u)
    vorticity = zeros(nx, ny)
    
    # Interior points using central differences
    for i in 2:nx-1
        for j in 2:ny-1
            dvdx = (result.v[i+1, j] - result.v[i-1, j]) / (2 * (result.x[2] - result.x[1]))
            dudy = (result.u[i, j+1] - result.u[i, j-1]) / (2 * (result.y[2] - result.y[1]))
            vorticity[i, j] = dvdx - dudy
        end
    end
    
    return vorticity
end

"""
    create_coordinate_meshgrid(result::NSResult) -> (Array, Array)

Create coordinate meshgrids for plotting.
"""
function create_coordinate_meshgrid(result::NSResult)
    X = result.x * ones(1, length(result.y))
    Y = ones(length(result.x), 1) * result.y'
    return X, Y
end

# Main visualization functions

"""
    plot_ns_solution(result::NSResult; backend=:plots, save_path=nothing, kwargs...)

Plot the complete Navier-Stokes solution including velocity and pressure fields.

# Arguments
- `result::NSResult`: Solution result from NSEMSolver
- `backend::Symbol`: Plotting backend (`:plots`, `:plotlyjs`, `:makie`, `:text`)
- `save_path::Union{String,Nothing}`: Optional path to save the plot
- `kwargs...`: Additional keyword arguments passed to plotting backend

# Returns
- Plot object (backend-dependent) or nothing for text output
"""
function plot_ns_solution(result::NSResult; backend=:plots, save_path=nothing, kwargs...)
    if backend == :text || (!HAS_PLOTS && !HAS_PLOTLYJS && !HAS_MAKIE)
        return _plot_text_summary(result)
    end
    
    if backend == :plots && HAS_PLOTS
        return _plot_with_plots(result, save_path; kwargs...)
    elseif backend == :plotlyjs && HAS_PLOTLYJS
        return _plot_with_plotlyjs(result, save_path; kwargs...)
    elseif backend == :makie && HAS_MAKIE
        return _plot_with_makie(result, save_path; kwargs...)
    else
        @warn "Requested backend $backend not available, falling back to text summary"
        return _plot_text_summary(result)
    end
end

"""
    plot_velocity_field(result::NSResult, component=:magnitude; backend=:plots, save_path=nothing, kwargs...)

Plot velocity field with options for different components.

# Arguments
- `result::NSResult`: Solution result
- `component::Symbol`: Component to plot (`:magnitude`, `:u`, `:v`, `:w`, `:streamlines`)
- `backend::Symbol`: Plotting backend
- `save_path::Union{String,Nothing}`: Optional path to save the plot
- `kwargs...`: Additional plotting arguments
"""
function plot_velocity_field(result::NSResult, component=:magnitude; backend=:plots, save_path=nothing, kwargs...)
    if backend == :text || !HAS_PLOTS
        println("Velocity Field Summary (component: $component):")
        if component == :magnitude
            vel_mag = compute_velocity_magnitude(result)
            println("  Magnitude: min=$(minimum(vel_mag)), max=$(maximum(vel_mag))")
        elseif component == :u
            println("  U-velocity: min=$(minimum(result.u)), max=$(maximum(result.u))")
        elseif component == :v
            println("  V-velocity: min=$(minimum(result.v)), max=$(maximum(result.v))")
        elseif component == :w && result.w !== nothing
            println("  W-velocity: min=$(minimum(result.w)), max=$(maximum(result.w))")
        end
        return nothing
    end
    
    p = _plot_velocity_component(result, component, backend; kwargs...)
    
    if save_path !== nothing && p !== nothing
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        Plots = PLOTS_MODULE
        Plots.savefig(p, joinpath(output_dir, save_path))
    end
    
    return p
end

"""
    plot_pressure_contours(result::NSResult; backend=:plots, levels=10, save_path=nothing, kwargs...)

Plot pressure field as contour plots.

# Arguments
- `result::NSResult`: Solution result
- `backend::Symbol`: Plotting backend
- `levels::Int`: Number of contour levels
- `save_path::Union{String,Nothing}`: Optional path to save the plot
- `kwargs...`: Additional plotting arguments
"""
function plot_pressure_contours(result::NSResult; backend=:plots, levels=10, save_path=nothing, kwargs...)
    if backend == :text || !HAS_PLOTS
        println("Pressure Field Summary:")
        println("  Pressure: min=$(minimum(result.p)), max=$(maximum(result.p))")
        println("  Requested contour levels: $levels")
        return nothing
    end
    
    p = _plot_pressure_field(result, backend, levels; kwargs...)
    
    if save_path !== nothing && p !== nothing
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        Plots = PLOTS_MODULE
        Plots.savefig(p, joinpath(output_dir, save_path))
    end
    
    return p
end

"""
    plot_streamlines(result::NSResult; backend=:plots, density=1.0, save_path=nothing, kwargs...)

Plot streamlines of the velocity field.

# Arguments
- `result::NSResult`: Solution result
- `backend::Symbol`: Plotting backend
- `density::Float64`: Streamline density
- `save_path::Union{String,Nothing}`: Optional path to save the plot
- `kwargs...`: Additional plotting arguments
"""
function plot_streamlines(result::NSResult; backend=:plots, density=1.0, save_path=nothing, kwargs...)
    if result.w !== nothing
        @warn "3D streamlines not implemented yet"
        return nothing
    end
    
    if backend == :text || !HAS_PLOTS
        println("Streamlines Summary:")
        println("  2D streamlines for velocity field")
        println("  Density: $density")
        return nothing
    end
    
    p = _plot_streamlines_2d(result, backend, density; kwargs...)
    
    if save_path !== nothing && p !== nothing
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        Plots = PLOTS_MODULE
        Plots.savefig(p, joinpath(output_dir, save_path))
    end
    
    return p
end

"""
    plot_convergence_history(result::NSResult; backend=:plots, log_scale=true, save_path=nothing, kwargs...)

Plot convergence history from the solution process.

# Arguments
- `result::NSResult`: Solution result
- `backend::Symbol`: Plotting backend
- `log_scale::Bool`: Use logarithmic scale for residuals
- `save_path::Union{String,Nothing}`: Optional path to save the plot
- `kwargs...`: Additional plotting arguments
"""
function plot_convergence_history(result::NSResult; backend=:plots, log_scale=true, save_path=nothing, kwargs...)
    if isempty(result.convergence_history)
        @warn "No convergence history available"
        return nothing
    end
    
    if backend == :text || !HAS_PLOTS
        println("Convergence History Summary:")
        println("  Total iterations: $(length(result.convergence_history))")
        println("  Final residual: $(@sprintf("%.2e", result.residual_norm))")
        println("  Converged: $(result.converged)")
        return nothing
    end
    
    p = _plot_convergence(result, backend, log_scale; kwargs...)
    
    if save_path !== nothing && p !== nothing
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        Plots = PLOTS_MODULE
        Plots.savefig(p, joinpath(output_dir, save_path))
    end
    
    return p
end

"""
    animate_solution(results::Vector{NSResult}, filename::String; backend=:plots, fps=10, kwargs...)

Create an animation from a time series of solutions.

# Arguments
- `results::Vector{NSResult}`: Time series of solution results
- `filename::String`: Output filename for animation
- `backend::Symbol`: Plotting backend
- `fps::Int`: Frames per second
- `kwargs...`: Additional plotting arguments
"""
function animate_solution(results::Vector{NSResult}, filename::String; backend=:plots, fps=10, kwargs...)
    if isempty(results)
        @warn "No results provided for animation"
        return nothing
    end
    
    if backend == :text || !HAS_PLOTS
        println("Animation Summary:")
        println("  Number of frames: $(length(results))")
        println("  Target FPS: $fps")
        println("  Output file: $filename")
        return nothing
    end
    
    return _create_animation(results, filename, backend, fps; kwargs...)
end

"""
    export_vtk(result::NSResult, filename::String; include_vorticity=true, kwargs...)

Export solution to VTK format for visualization in ParaView or other tools.

# Arguments
- `result::NSResult`: Solution result
- `filename::String`: Output filename (without extension)
- `include_vorticity::Bool`: Include vorticity field in export
- `kwargs...`: Additional export options
"""
function export_vtk(result::NSResult, filename::String; include_vorticity=true, kwargs...)
    if !HAS_WRITEVTK
        @warn "WriteVTK.jl not available. Cannot export VTK files."
        return _export_csv(result, filename * ".csv")
    end
    
    return _export_vtk_impl(result, filename, include_vorticity; kwargs...)
end

# Backend-specific implementation functions

function _plot_text_summary(result::NSResult)
    println("Navier-Stokes Solution Summary:")
    println("  Grid size: $(size(result.u))")
    println("  Velocity u: min=$(@sprintf("%.3f", minimum(result.u))), max=$(@sprintf("%.3f", maximum(result.u)))")
    println("  Velocity v: min=$(@sprintf("%.3f", minimum(result.v))), max=$(@sprintf("%.3f", maximum(result.v)))")
    if result.w !== nothing
        println("  Velocity w: min=$(@sprintf("%.3f", minimum(result.w))), max=$(@sprintf("%.3f", maximum(result.w)))")
    end
    println("  Pressure p: min=$(@sprintf("%.3f", minimum(result.p))), max=$(@sprintf("%.3f", maximum(result.p)))")
    
    vel_mag = compute_velocity_magnitude(result)
    println("  Velocity magnitude: min=$(@sprintf("%.3f", minimum(vel_mag))), max=$(@sprintf("%.3f", maximum(vel_mag)))")
    
    println("  Converged: $(result.converged)")
    println("  Iterations: $(result.iterations)")
    println("  Final residual: $(@sprintf("%.2e", result.residual_norm))")
    println("  Solve time: $(@sprintf("%.3f", result.solve_time)) seconds")
    
    return nothing
end

# Backend-specific plotting implementations

function _plot_with_plots(result::NSResult, save_path; kwargs...)
    if !HAS_PLOTS
        return _plot_text_summary(result)
    end
    
    Plots = PLOTS_MODULE
    
    # Create subplots for velocity components and pressure
    X, Y = create_coordinate_meshgrid(result)
    vel_mag = compute_velocity_magnitude(result)
    
    # Create 2x2 subplot layout
    p1 = Plots.heatmap(result.x, result.y, result.u', 
                       title="U-velocity", xlabel="x", ylabel="y", 
                       aspect_ratio=:equal, color=:viridis)
    
    p2 = Plots.heatmap(result.x, result.y, result.v', 
                       title="V-velocity", xlabel="x", ylabel="y", 
                       aspect_ratio=:equal, color=:viridis)
    
    p3 = Plots.heatmap(result.x, result.y, vel_mag', 
                       title="Velocity Magnitude", xlabel="x", ylabel="y", 
                       aspect_ratio=:equal, color=:plasma)
    
    p4 = Plots.contour(result.x, result.y, result.p', 
                       title="Pressure Contours", xlabel="x", ylabel="y", 
                       aspect_ratio=:equal, levels=10)
    
    plot_combined = Plots.plot(p1, p2, p3, p4, layout=(2,2), size=(800, 600))
    
    if save_path !== nothing
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        Plots.savefig(plot_combined, joinpath(output_dir, save_path))
    end
    
    return plot_combined
end

function _plot_with_plotlyjs(result::NSResult, save_path; kwargs...)
    if !HAS_PLOTLYJS
        return _plot_text_summary(result)
    end
    
    # Implementation for PlotlyJS backend
    # This would create interactive plots
    @warn "PlotlyJS backend implementation in progress"
    return _plot_text_summary(result)
end

function _plot_with_makie(result::NSResult, save_path; kwargs...)
    if !HAS_MAKIE
        return _plot_text_summary(result)
    end
    
    # Implementation for Makie backend
    # This would create high-quality 3D plots
    @warn "Makie backend implementation in progress"
    return _plot_text_summary(result)
end

function _plot_velocity_component(result::NSResult, component::Symbol, backend::Symbol; kwargs...)
    if !HAS_PLOTS && backend == :plots
        return _plot_text_summary(result)
    end
    
    Plots = PLOTS_MODULE
    X, Y = create_coordinate_meshgrid(result)
    
    if component == :magnitude
        data = compute_velocity_magnitude(result)'
        title = "Velocity Magnitude"
        color = :plasma
    elseif component == :u
        data = result.u'
        title = "U-velocity Component"
        color = :viridis
    elseif component == :v
        data = result.v'
        title = "V-velocity Component"
        color = :viridis
    elseif component == :w && result.w !== nothing
        data = result.w'
        title = "W-velocity Component"
        color = :viridis
    else
        @warn "Unknown velocity component: $component"
        return nothing
    end
    
    p = Plots.heatmap(result.x, result.y, data,
                      title=title, xlabel="x", ylabel="y",
                      aspect_ratio=:equal, color=color)
    
    if component == :streamlines
        # Add streamlines overlay
        skip = max(1, length(result.x) ÷ 20)  # Downsample for streamlines
        x_stream = result.x[1:skip:end]
        y_stream = result.y[1:skip:end]
        u_stream = result.u[1:skip:end, 1:skip:end]
        v_stream = result.v[1:skip:end, 1:skip:end]
        
        # Simple streamline approximation using quiver plot
        Plots.quiver!(x_stream, y_stream, quiver=(u_stream, v_stream), 
                      color=:white, alpha=0.7)
    end
    
    return p
end

function _plot_pressure_field(result::NSResult, backend::Symbol, levels::Int; kwargs...)
    if !HAS_PLOTS && backend == :plots
        return _plot_text_summary(result)
    end
    
    Plots = PLOTS_MODULE
    
    p = Plots.contourf(result.x, result.y, result.p',
                       title="Pressure Field", xlabel="x", ylabel="y",
                       aspect_ratio=:equal, levels=levels, color=:RdYlBu)
    
    return p
end

function _plot_streamlines_2d(result::NSResult, backend::Symbol, density::Float64; kwargs...)
    if !HAS_PLOTS && backend == :plots
        return _plot_text_summary(result)
    end
    
    Plots = PLOTS_MODULE
    
    # Create background velocity magnitude
    vel_mag = compute_velocity_magnitude(result)
    p = Plots.heatmap(result.x, result.y, vel_mag',
                      title="Streamlines on Velocity Magnitude", 
                      xlabel="x", ylabel="y",
                      aspect_ratio=:equal, color=:plasma, alpha=0.7)
    
    # Add streamlines (simplified implementation)
    skip = max(1, round(Int, length(result.x) / (density * 20)))
    x_stream = result.x[1:skip:end]
    y_stream = result.y[1:skip:end]
    u_stream = result.u[1:skip:end, 1:skip:end]
    v_stream = result.v[1:skip:end, 1:skip:end]
    
    Plots.quiver!(x_stream, y_stream, quiver=(u_stream * 0.1, v_stream * 0.1), 
                  color=:white, alpha=0.8, linewidth=1)
    
    return p
end

function _plot_convergence(result::NSResult, backend::Symbol, log_scale::Bool; kwargs...)
    if !HAS_PLOTS && backend == :plots
        return _plot_text_summary(result)
    end
    
    Plots = PLOTS_MODULE
    
    iterations = 1:length(result.convergence_history)
    
    if log_scale
        p = Plots.plot(iterations, result.convergence_history,
                       title="Convergence History", xlabel="Iteration", ylabel="Residual",
                       yscale=:log10, linewidth=2, marker=:circle, markersize=3,
                       label="Residual", legend=:topright)
    else
        p = Plots.plot(iterations, result.convergence_history,
                       title="Convergence History", xlabel="Iteration", ylabel="Residual",
                       linewidth=2, marker=:circle, markersize=3,
                       label="Residual", legend=:topright)
    end
    
    # Add tolerance line if available
    if hasfield(typeof(result.options), :tol)
        Plots.hline!([result.options.tol], linestyle=:dash, 
                     color=:red, label="Tolerance")
    end
    
    return p
end

function _create_animation(results::Vector{NSResult}, filename::String, backend::Symbol, fps::Int; kwargs...)
    if !HAS_PLOTS && backend == :plots
        @warn "Animation requires Plots.jl backend"
        return nothing
    end
    
    Plots = PLOTS_MODULE
    
    # Create animation of velocity magnitude evolution
    anim = Plots.Animation()
    
    for (i, result) in enumerate(results)
        vel_mag = compute_velocity_magnitude(result)
        
        p = Plots.heatmap(result.x, result.y, vel_mag',
                          title="Velocity Evolution - Frame $i",
                          xlabel="x", ylabel="y",
                          aspect_ratio=:equal, color=:plasma,
                          clims=(0, maximum([maximum(compute_velocity_magnitude(r)) for r in results])))
        
        Plots.frame(anim, p)
    end
    
    # Save animation
    output_dir = "../visualization_outputs"
    mkpath(output_dir)
    output_path = joinpath(output_dir, filename)
    Plots.gif(anim, output_path, fps=fps)
    @info "Animation saved to $output_path"
    
    return anim
end

# VTK export implementation

function _export_vtk_impl(result::NSResult, filename::String, include_vorticity::Bool; kwargs...)
    if !HAS_WRITEVTK
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        return _export_csv(result, joinpath(output_dir, filename * ".csv"))
    end
    
    WriteVTK = WRITEVTK_MODULE
    
    # Prepare output path
    output_dir = "../visualization_outputs"
    mkpath(output_dir)
    output_filename = joinpath(output_dir, filename)
    
    # Create structured grid
    nx, ny = length(result.x), length(result.y)
    
    if result.w !== nothing
        # 3D case
        nz = length(result.z)
        vtk_grid = WriteVTK.vtk_grid(output_filename, result.x, result.y, result.z)
    else
        # 2D case - create thin 3D grid
        z_coords = [0.0, 0.1]  # Thin layer for 2D visualization
        vtk_grid = WriteVTK.vtk_grid(output_filename, result.x, result.y, z_coords)
    end
    
    # Add velocity field
    if result.w !== nothing
        # 3D velocity
        velocity = zeros(3, nx, ny, length(result.z))
        velocity[1, :, :, :] .= result.u
        velocity[2, :, :, :] .= result.v
        velocity[3, :, :, :] .= result.w
    else
        # 2D velocity (extended to thin 3D)
        velocity = zeros(3, nx, ny, 2)
        for k in 1:2
            velocity[1, :, :, k] .= result.u
            velocity[2, :, :, k] .= result.v
            velocity[3, :, :, k] .= 0.0
        end
    end
    
    vtk_grid["velocity"] = velocity
    
    # Add pressure field
    if result.w !== nothing
        vtk_grid["pressure"] = result.p
    else
        pressure_3d = zeros(nx, ny, 2)
        for k in 1:2
            pressure_3d[:, :, k] .= result.p
        end
        vtk_grid["pressure"] = pressure_3d
    end
    
    # Add velocity magnitude
    vel_mag = compute_velocity_magnitude(result)
    if result.w === nothing
        vel_mag_3d = zeros(nx, ny, 2)
        for k in 1:2
            vel_mag_3d[:, :, k] .= vel_mag
        end
        vtk_grid["velocity_magnitude"] = vel_mag_3d
    else
        vtk_grid["velocity_magnitude"] = vel_mag
    end
    
    # Add vorticity if requested
    if include_vorticity && result.w === nothing  # 2D case only for now
        vort = compute_vorticity(result)
        vort_3d = zeros(nx, ny, 2)
        for k in 1:2
            vort_3d[:, :, k] .= vort
        end
        vtk_grid["vorticity"] = vort_3d
    end
    
    # Save the file
    WriteVTK.vtk_save(vtk_grid)
    @info "VTK file saved as $(output_filename).vtu"
    
    return vtk_grid
end

function _export_csv(result::NSResult, filename::String)
    @info "Exporting to CSV format as fallback"
    
    # Handle output directory if filename doesn't already include it
    if !startswith(filename, "../visualization_outputs")
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        filename = joinpath(output_dir, filename)
    end
    
    try
        open(filename, "w") do io
            if result.w !== nothing
                # 3D case
                println(io, "x,y,z,u,v,w,p,vel_magnitude")
                for k in 1:length(result.z), j in 1:length(result.y), i in 1:length(result.x)
                    if i <= size(result.u, 1) && j <= size(result.u, 2) && k <= size(result.u, 3)
                        vel_mag = sqrt(result.u[i,j,k]^2 + result.v[i,j,k]^2 + result.w[i,j,k]^2)
                        println(io, "$(result.x[i]),$(result.y[j]),$(result.z[k]),$(result.u[i,j,k]),$(result.v[i,j,k]),$(result.w[i,j,k]),$(result.p[i,j,k]),$vel_mag")
                    end
                end
            else
                # 2D case
                println(io, "x,y,u,v,p,vel_magnitude")
                for j in 1:length(result.y), i in 1:length(result.x)
                    if i <= size(result.u, 1) && j <= size(result.u, 2)
                        vel_mag = sqrt(result.u[i,j]^2 + result.v[i,j]^2)
                        println(io, "$(result.x[i]),$(result.y[j]),$(result.u[i,j]),$(result.v[i,j]),$(result.p[i,j]),$vel_mag")
                    end
                end
            end
        end
        @info "Solution exported to $filename"
    catch e
        @error "Failed to export solution: $e"
    end
    
    return nothing
end

# Utility functions for specialized visualizations

"""
    plot_domain_refinement(result::NSResult; backend=:plots, save_path=nothing, kwargs...)

Visualize adaptive refinement levels across the computational domain.
"""
function plot_domain_refinement(result::NSResult; backend=:plots, save_path=nothing, kwargs...)
    if result.multidomain === nothing || 
       !hasfield(typeof(result.multidomain), :refinement_map) ||
       result.multidomain.refinement_map === nothing
        @warn "No refinement information available"
        return nothing
    end
    
    if backend == :text || !HAS_PLOTS
        println("Domain Refinement Summary:")
        println("  Refinement map size: $(size(result.multidomain.refinement_map))")
        println("  Refinement levels: $(minimum(result.multidomain.refinement_map)) to $(maximum(result.multidomain.refinement_map))")
        return nothing
    end
    
    Plots = PLOTS_MODULE
    
    p = Plots.heatmap(result.multidomain.refinement_map,
                      title="Adaptive Refinement Levels",
                      xlabel="Block X", ylabel="Block Y",
                      aspect_ratio=:equal, color=:viridis)
    
    if save_path !== nothing && p !== nothing
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        Plots.savefig(p, joinpath(output_dir, save_path))
    end
    
    return p
end

"""
    plot_velocity_vectors(result::NSResult; backend=:plots, skip=5, scale=1.0, save_path=nothing, kwargs...)

Plot velocity field as vector arrows.

# Arguments
- `result::NSResult`: Solution result
- `backend::Symbol`: Plotting backend
- `skip::Int`: Skip factor for vector density (plot every skip-th vector)
- `scale::Float64`: Arrow scaling factor
- `save_path::Union{String,Nothing}`: Optional path to save the plot
- `kwargs...`: Additional plotting arguments
"""
function plot_velocity_vectors(result::NSResult; backend=:plots, skip=5, scale=1.0, save_path=nothing, kwargs...)
    if backend == :text || !HAS_PLOTS
        println("Velocity Vectors Summary:")
        println("  Vector field for velocity components u, v")
        println("  Skip factor: $skip")
        println("  Scale factor: $scale")
        return nothing
    end
    
    Plots = PLOTS_MODULE
    
    # Create background with velocity magnitude
    vel_mag = compute_velocity_magnitude(result)
    p = Plots.heatmap(result.x, result.y, vel_mag',
                      title="Velocity Vector Field",
                      xlabel="x", ylabel="y",
                      aspect_ratio=:equal, color=:plasma, alpha=0.6)
    
    # Add velocity vectors
    x_vec = result.x[1:skip:end]
    y_vec = result.y[1:skip:end]
    u_vec = result.u[1:skip:end, 1:skip:end] * scale
    v_vec = result.v[1:skip:end, 1:skip:end] * scale
    
    Plots.quiver!(x_vec, y_vec, quiver=(u_vec, v_vec),
                  color=:white, alpha=0.8, linewidth=1.5)
    
    if save_path !== nothing && p !== nothing
        output_dir = "../visualization_outputs"
        mkpath(output_dir)
        Plots.savefig(p, joinpath(output_dir, save_path))
    end
    
    return p
end

"""
    create_visualization_summary(result::NSResult) -> Dict

Create a summary dictionary of visualization data for external processing.
"""
function create_visualization_summary(result::NSResult)
    summary = Dict{String, Any}()
    
    # Basic solution info
    summary["grid_size"] = size(result.u)
    summary["is_3d"] = result.w !== nothing
    summary["converged"] = result.converged
    summary["iterations"] = result.iterations
    summary["residual_norm"] = result.residual_norm
    summary["solve_time"] = result.solve_time
    
    # Field statistics
    summary["velocity_stats"] = Dict(
        "u_min" => minimum(result.u),
        "u_max" => maximum(result.u),
        "v_min" => minimum(result.v),
        "v_max" => maximum(result.v)
    )
    
    if result.w !== nothing
        summary["velocity_stats"]["w_min"] = minimum(result.w)
        summary["velocity_stats"]["w_max"] = maximum(result.w)
    end
    
    summary["pressure_stats"] = Dict(
        "p_min" => minimum(result.p),
        "p_max" => maximum(result.p)
    )
    
    # Derived quantities
    vel_mag = compute_velocity_magnitude(result)
    summary["velocity_magnitude_stats"] = Dict(
        "min" => minimum(vel_mag),
        "max" => maximum(vel_mag),
        "mean" => sum(vel_mag) / length(vel_mag)
    )
    
    # Coordinate ranges
    summary["coordinate_ranges"] = Dict(
        "x_range" => [minimum(result.x), maximum(result.x)],
        "y_range" => [minimum(result.y), maximum(result.y)]
    )
    
    if result.z !== nothing
        summary["coordinate_ranges"]["z_range"] = [minimum(result.z), maximum(result.z)]
    end
    
    # Convergence info
    if !isempty(result.convergence_history)
        summary["convergence_info"] = Dict(
            "initial_residual" => result.convergence_history[1],
            "final_residual" => result.convergence_history[end],
            "reduction_factor" => result.convergence_history[1] / result.convergence_history[end]
        )
    end
    
    return summary
end

# Backwards compatibility aliases
const export_solution = export_vtk  # Maintain compatibility with old function name